/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package genomemap.worker.helper;

import javautil.lang.MathUtil;

/**
 * Class description
 *
 * @version        Enter version here..., 12/11/23
 * @author         Susanta Tewari
 */
public class StdErrHelper {

    public static double[][] getInfoDerivs(int[][] genotype, int[] genotypeFreq,
            double[] cProbFinal, int distinctGenotypes) {

        int loci           = cProbFinal.length + 1;
        int height         = 3;
        boolean verbose    = false;
        double[] probFOld1 = new double[distinctGenotypes];

        /**
         * CREATING activeKFirst,cArrayFirst,linkInfoFirst,kArrayFirst
         */
        MyInteger activeKFirstIntObj = new MyInteger(0);

        GHoodHelper.getActiveXoverFirst(height, genotype, activeKFirstIntObj);

        int activeKFirst      = activeKFirstIntObj.value;
        int[][] kArrayFirst   = new int[activeKFirst][height - 1];
        int[] linkInfoFirst   = new int[activeKFirst];
        int[][][] cArrayFirst = new int[activeKFirst][distinctGenotypes][4];

        GHoodHelper.getInfo4First(height, genotype, kArrayFirst, linkInfoFirst, cArrayFirst);

        /**
         * CREATING activeKLast,cArrayLast,kArrayLast
         */
        int[] activeKLast = new int[6];

        GHoodHelper.getActiveXoverLastArray(height, genotype, activeKLast);

        int[][][] kArrayLast   = new int[6][][];
        int[][][][] cArrayLast = new int[6][][][];

        for (int i = 0; i < 6; i++) {

            kArrayLast[i] = new int[activeKLast[i]][height - 1];
            cArrayLast[i] = new int[activeKLast[i]][distinctGenotypes][4];
        }

        GHoodHelper.getInfo4Last(height, genotype, kArrayLast, cArrayLast);

        /**
         * Create the single and double derivatives for the information matrix
         */

        /**
         * FIRST PART:
         */
        double[] lastCProb = GHoodHelper.ChopAtoB(cProbFinal, loci - height + 1, loci - 1);
        int[][] spores                  = GHoodHelper.getSpores();
        double probDeri                 = 0.0,
               probSingle               = 0.0,
               probSingle2              = 0.0;
        double[][][][][] FirstTempInfo =
            new double[6][11][height - 1][height - 1][distinctGenotypes];
        double[][][][][] FirstTempInfo2 =
            new double[6][11][height - 1][height - 1][distinctGenotypes];
        double[][][][] FirstTempSingleD = new double[6][11][height - 1][distinctGenotypes];
        double[][][] FirstTempSum       = new double[6][distinctGenotypes][11];

        for (int index1 = 0; index1 < 6; index1++) {

            for (int index6 = 0; index6 < activeKLast[index1]; index6++) {

                double prob = GHoodHelper.kProb(lastCProb, kArrayLast[index1][index6]);

                for (int index2 = 0; index2 < 11; index2++) {

                    for (int index5 = 0; index5 < distinctGenotypes; index5++) {

                        double condProb = MathUtil.Sum(
                                              GHoodHelper.UpdateSpore(
                                                  cArrayLast[index1][index6][index5],
                                                  spores[index2])) / 4.0;

                        for (int index3 = loci - height; index3 < loci - 1; index3++) {

                            probSingle = kProbSingle(lastCProb, kArrayLast[index1][index6],
                                                     index3 - loci + height + 1);

                            for (int index4 = loci - height; index4 < loci - 1; index4++) {

                                if (index3 != index4) {

                                    probDeri = kProbDeriOff(lastCProb, kArrayLast[index1][index6],
                                                            index3 - loci + height + 1,
                                                            index4 - loci + height + 1);

                                } else {

                                    probDeri = kProbDeriDiag(lastCProb, kArrayLast[index1][index6],
                                                             index3 - loci + height + 1);
                                }

                                FirstTempInfo[index1][index2][index3 - loci + height][index4 - loci + height][index5] +=
                                    condProb * probDeri;
                                probSingle2 = kProbSingle(lastCProb, kArrayLast[index1][index6],
                                                          index4 - loci + height + 1);
                                FirstTempInfo2[index1][index2][index3 - loci + height][index4 - loci + height][index5] +=
                                    condProb / prob * probSingle * probSingle2;
                            }

                            FirstTempSingleD[index1][index2][index3 - loci + height][index5] +=
                                condProb * probSingle;


                            // FirstTempSingleD2[index1][index2][index3-loci+height][index5]+=condProb*probSingle;
                        }

                        FirstTempSum[index1][index5][index2] += condProb * prob;
                    }
                }
            }
        }

        /**
         * RECURSIVELY FOLLOW ALL THE SEGMENTS
         */
        int startRowIndex           = 0;
        int endRowIndex = loci - height;    // note that there is a common gene in each interval
        int startPos                = 0;
        int endPos                  = loci - height;
        int segments                = GHoodHelper.getSegments(loci, height);
        double[][][][][] temp1Info  = FirstTempInfo.clone();
        double[][][][][] temp1Info2 = FirstTempInfo2.clone();
        double[][][][] temp1SingleD = FirstTempSingleD.clone();
        double[][][] temp1Sum       = FirstTempSum.clone();

        FirstTempInfo    = null;
        FirstTempInfo2   = null;
        FirstTempSingleD = null;
        FirstTempSum     = null;

        for (int indexS = 0; indexS < segments - 2; indexS++) {

            startPos = endPos - height + 2;

            double[][][][][] temp2Info =
                new double[6][11][loci - startPos][loci - startPos][distinctGenotypes];
            double[][][][][] temp2Info2 =
                new double[6][11][loci - startPos][loci - startPos][distinctGenotypes];
            double[][][][] temp2SingleD = new double[6][11][loci - startPos][distinctGenotypes];
            double[][][] temp2Sum       = new double[6][distinctGenotypes][11];
            double[] tempCProb          = GHoodHelper.ChopAtoB(cProbFinal, startPos, endPos);


            // setting up the sth temporary segment
            int[] activeKTemp = new int[6];


            // Get the relevant data
            startRowIndex = endRowIndex - height + 1;    // endRowIndex is changed in the bottom.

            int[][] data = GHoodHelper.getGenotypeSegment(startRowIndex, endRowIndex, genotype);


            // Calculate active crossovers. Stored in activeKTemp.
            GHoodHelper.ActiveKLoopTemp(0, height, indexS, activeKTemp, data);    // note this uses data.


            // allocate cells.
            int[][][][] cArrayTemp = new int[6][][][];
            int[][] linkInfoTemp   = new int[6][];
            int[][][] kArrayTemp   = new int[6][][];

            for (int i = 0; i < 6; i++) {

                cArrayTemp[i]   = new int[activeKTemp[i]][distinctGenotypes][4];
                linkInfoTemp[i] = new int[activeKTemp[i]];
                kArrayTemp[i]   = new int[activeKTemp[i]][height - 1];
            }


            // run the loop to set kArrayTemp,cArrayTemp,linkInfoTemp.
            GHoodHelper.LoopTemp(0, height, indexS, data, kArrayTemp, linkInfoTemp, cArrayTemp);    // note this uses data


            // linking at work
            for (int index1 = 0; index1 < 6; index1++) {

                for (int indexK = 0; indexK < activeKTemp[index1]; indexK++) {

                    double prob = GHoodHelper.kProb(tempCProb, kArrayTemp[index1][indexK]);

                    for (int index2 = 0; index2 < 11; index2++) {

                        for (int index5 = 0; index5 < distinctGenotypes; index5++) {

                            int spike = GHoodHelper.getSporeMatch(
                                            GHoodHelper.UpdateSpore(
                                                cArrayTemp[index1][indexK][index5],
                                                spores[index2]));
                            int index01 = linkInfoTemp[index1][indexK];

                            temp2Sum[index1][index5][index2] += temp1Sum[index01][index5][spike]
                                    * prob;

                            for (int index3 = startPos - 1; index3 < loci - 1; index3++) {

                                for (int index4 = startPos - 1; index4 < loci - 1; index4++) {

                                    if (((index3 >= startPos - 1) && (index3 <= endPos - 1))
                                            && ((index4 >= startPos - 1)
                                                && (index4 <= endPos - 1))) {

                                        if (index3 != index4) {

                                            probDeri = kProbDeriOff(tempCProb,
                                                                    kArrayTemp[index1][indexK],
                                                                    index3 - startPos + 1 + 1,
                                                                    index4 - startPos + 1 + 1);

                                        } else {

                                            probDeri = kProbDeriDiag(tempCProb,
                                                                     kArrayTemp[index1][indexK],
                                                                     index3 - startPos + 1 + 1);
                                        }

                                        temp2Info[index1][index2][index3 - startPos + 1][index4 - startPos + 1][index5] +=
                                            temp1Sum[index01][index5][spike] * probDeri;
                                        probSingle = kProbSingle(tempCProb,
                                                                 kArrayTemp[index1][indexK],
                                                                 index3 - startPos + 1 + 1);
                                        probSingle2 = kProbSingle(tempCProb,
                                                                  kArrayTemp[index1][indexK],
                                                                  index4 - startPos + 1 + 1);
                                        temp2Info2[index1][index2][index3 - startPos + 1][index4 - startPos + 1][index5] +=
                                            temp1Sum[index01][index5][spike] * probSingle
                                            * probSingle2 / prob;

                                    } else {

                                        if (((index3 >= startPos - 1) && (index3 <= endPos - 1))
                                                && ((index4 >= endPos) && (index4 <= loci - 2))) {

                                            probSingle = kProbSingle(tempCProb,
                                                                     kArrayTemp[index1][indexK],
                                                                     index3 - startPos + 1 + 1);
                                            temp2Info[index1][index2][index3 - startPos + 1][index4 - startPos + 1][index5] +=
                                                temp1SingleD[index01][spike][index4 - endPos][index5]
                                                * probSingle;
                                            temp2Info2[index1][index2][index3 - startPos + 1][index4 - startPos + 1][index5] +=
                                                temp1SingleD[index01][spike][index4 - endPos][index5]
                                                * probSingle;

                                        } else {

                                            if (((index3 >= endPos) && (index3 <= loci - 2))
                                                    && ((index4 >= startPos - 1)
                                                        && (index3 <= endPos - 1))) {

                                                probSingle =
                                                    kProbSingle(tempCProb,
                                                                kArrayTemp[index1][indexK],
                                                                index4 - startPos + 1 + 1);
                                                temp2Info[index1][index2][index3 - startPos + 1][index4 - startPos + 1][index5] +=
                                                    temp1SingleD[index01][spike][index3 - endPos][index5]
                                                    * probSingle;
                                                temp2Info2[index1][index2][index3 - startPos + 1][index4 - startPos + 1][index5] +=
                                                    temp1SingleD[index01][spike][index3 - endPos][index5]
                                                    * probSingle;

                                            } else {

                                                if (((index3 >= endPos) && (index3 <= loci - 2))
                                                        && ((index4 >= endPos)
                                                            && (index4 <= loci - 2))) {

                                                    temp2Info[index1][index2][index3 - startPos + 1][index4 - startPos + 1][index5] +=
                                                        temp1Info[index01][spike][index3 - endPos][index4 - endPos][index5]
                                                        * prob;
                                                    temp2Info2[index1][index2][index3 - startPos + 1][index4 - startPos + 1][index5] +=
                                                        temp1Info2[index01][spike][index3 - endPos][index4 - endPos][index5]
                                                        * prob;
                                                }
                                            }
                                        }
                                    }

                                    if ((index3 >= startPos - 1) && (index3 <= endPos - 1)) {

                                        probSingle = kProbSingle(tempCProb,
                                                                 kArrayTemp[index1][indexK],
                                                                 index3 - startPos + 1 + 1);
                                        temp2SingleD[index1][index2][index3 - startPos + 1][index5] +=
                                            temp1Sum[index01][index5][spike] * probSingle;

                                    } else {

                                        temp2SingleD[index1][index2][index3 - startPos + 1][index5] +=
                                            temp1SingleD[index01][spike][index3 - endPos][index5]
                                            * prob;
                                    }
                                }
                            }
                        }
                    }
                }
            }


            // change temp1 to temp2
            temp1Info    = temp2Info.clone();
            temp1Info2   = temp2Info2.clone();
            temp1SingleD = temp2SingleD.clone();
            temp1Sum     = temp2Sum.clone();
            endPos       = startPos - 1;
            endRowIndex  = startRowIndex;

            if (verbose) System.out.println("From StdErrHelper: Segment # passed " + indexS);
        }

        /*
         * ---------------------------------------------
         * The Final Loop.
         * ---------------------------------------------
         */
        double[] firstCProb    = GHoodHelper.ChopAtoB(cProbFinal, 1, height - 1);
        double[][][] tempInfo  = new double[loci - 1][loci - 1][distinctGenotypes];
        double[][][] tempInfo2 = new double[loci - 1][loci - 1][distinctGenotypes];
        double[][] tempSingleD = new double[loci - 1][distinctGenotypes];

        for (int index1 = 0; index1 < activeKFirst; index1++) {

            int[] spike = new int[distinctGenotypes];

            for (int i = 0; i < distinctGenotypes; i++) {

                spike[i] = GHoodHelper.getSporeMatch(cArrayFirst[index1][i]);
            }

            double prob = GHoodHelper.kProb(firstCProb, kArrayFirst[index1]);
            int index01 = linkInfoFirst[index1];

            for (int index2 = 0; index2 < distinctGenotypes; index2++) {

                probFOld1[index2] += temp1Sum[index01][index2][spike[index2]] * prob;

                for (int index3 = 0; index3 < loci - 1; index3++) {

                    for (int index4 = 0; index4 < loci - 1; index4++) {

                        if (((index3 >= 0) && (index3 <= height - 2))
                                && ((index4 >= 0) && (index4 <= height - 2))) {

                            if (index3 != index4) {

                                probDeri = kProbDeriOff(firstCProb, kArrayFirst[index1],
                                                        index3 + 1, index4 + 1);

                            } else {

                                probDeri = kProbDeriDiag(firstCProb, kArrayFirst[index1],
                                                         index3 + 1);
                            }

                            tempInfo[index3][index4][index2] +=
                                temp1Sum[index01][index2][spike[index2]] * probDeri;
                            probSingle = kProbSingle(firstCProb, kArrayFirst[index1], index3 + 1);
                            probSingle2 = kProbSingle(firstCProb, kArrayFirst[index1], index4 + 1);
                            tempInfo2[index3][index4][index2] +=
                                temp1Sum[index01][index2][spike[index2]] * probSingle * probSingle2
                                / prob;

                        } else {

                            if (((index3 >= 0) && (index3 <= height - 2))
                                    && ((index4 >= height - 1) && (index4 <= loci - 2))) {

                                probSingle = kProbSingle(firstCProb, kArrayFirst[index1],
                                                         index3 + 1);
                                tempInfo[index3][index4][index2] +=
                                    temp1SingleD[index01][spike[index2]][index4 - endPos][index2]
                                    * probSingle;
                                tempInfo2[index3][index4][index2] +=
                                    temp1SingleD[index01][spike[index2]][index4 - endPos][index2]
                                    * probSingle;

                            } else {

                                if (((index3 >= height - 1) && (index3 <= loci - 2))
                                        && ((index4 >= 0) && (index4 <= height - 2))) {

                                    probSingle = kProbSingle(firstCProb, kArrayFirst[index1],
                                                             index4 + 1);
                                    tempInfo[index3][index4][index2] +=
                                        temp1SingleD[index01][spike[index2]][index3 - endPos][index2]
                                        * probSingle;
                                    tempInfo2[index3][index4][index2] +=
                                        temp1SingleD[index01][spike[index2]][index3 - endPos][index2]
                                        * probSingle;

                                } else {

                                    if (((index3 >= height - 1) && (index3 <= loci - 2))
                                            && ((index4 >= height - 1) && (index4 <= loci - 2))) {

                                        tempInfo[index3][index4][index2] +=
                                            temp1Info[index01][spike[index2]][index3 - endPos][index4 - endPos][index2]
                                            * prob;
                                        tempInfo2[index3][index4][index2] +=
                                            temp1Info2[index01][spike[index2]][index3 - endPos][index4 - endPos][index2]
                                            * prob;
                                    }
                                }
                            }
                        }
                    }

                    if ((index3 >= 0) && (index3 <= height - 2)) {

                        probSingle = kProbSingle(firstCProb, kArrayFirst[index1], index3 + 1);
                        tempSingleD[index3][index2] += temp1Sum[index01][index2][spike[index2]]
                                                       * probSingle;

                    } else {

                        tempSingleD[index3][index2] +=
                            temp1SingleD[index01][spike[index2]][index3 - endPos][index2] * prob;
                    }
                }
            }
        }

        if (verbose) System.out.println("From StdErrHelper : Finshed all looping.  ");


        // Getting the Expected information matrix
        double[][] info  = new double[loci - 1][loci - 1];
        double[][] info1 = new double[loci - 1][loci - 1];
        double sum$cov   = 0.0;

        for (int index1 = 0; index1 < loci - 1; index1++) {

            for (int index2 = 0; index2 < loci - 1; index2++) {

                sum$cov = 0.0;

                for (int index3 = 0; index3 < distinctGenotypes; index3++) {


                    // sum$cov+=(probFOld1[index3]*tempInfo[index1][index2][index3]-tempSingleD[index1][index3]*tempSingleD[index2][index3])/probFOld1[index3];
                    sum$cov += (tempInfo2[index1][index2][index3]
                                - tempInfo[index1][index2][index3]) * genotypeFreq[index3]
                                    / probFOld1[index3];
                }


                // info[index1][index2]=sum$cov*(-1.0)*total_obs;
                info[index1][index2] = sum$cov;
            }
        }

        return info;
    }

    /**
     * Method description
     *
     * @param cProbOld description
     * @param k description
     * @param i description
     * @param j description
     *
     * @return description
     */
    private static double kProbDeriOff(double[] cProbOld, int[] k, int i, int j) {

        double product = 1;

        for (int index = 0; index < k.length; index++) {

            if ((index != (i - 1)) && (index != (j - 1))) {
                product *= GHoodHelper.lociProb(k[index], cProbOld[index]);
            } else {

                if (index == i - 1) {
                    product = product * probD(cProbOld, k, i);
                } else {
                    product = product * probD(cProbOld, k, j);
                }
            }
        }

        return product;
    }

    /**
     * Method description
     *
     * @param cProbOld description
     * @param k description
     * @param i description
     *
     * @return description
     */
    private static double probD(double[] cProbOld, int[] k, int i) {

        double result = 0.0;

        if (k[i - 1] == 0) {
            result = 2 * (1 - cProbOld[i - 1]) * (-1);
        } else if (((k[i - 1] >= 1) && (k[i - 1] <= 5)) || (k[i - 1] == 10) || (k[i - 1] == 15)
                   || (k[i - 1] == 20)) {
            result = 0.25 - cProbOld[i - 1] / 2.0;
        } else {
            result = cProbOld[i - 1] / 8.0;
        }

        return result;
    }

    /**
     * Method description
     *
     * @param cProbOld description
     * @param k description
     * @param i description
     *
     * @return description
     */
    private static double kProbDeriDiag(double[] cProbOld, int[] k, int i) {

        double product = 1;

        for (int index = 0; index < k.length; index++) {

            if (index != (i - 1)) {
                product *= GHoodHelper.lociProb(k[index], cProbOld[index]);
            } else {
                product = product * probD2(cProbOld, k, i);
            }
        }

        return product;
    }

    /**
     * Method description
     *
     * @param cProbOld description
     * @param k description
     * @param i description
     *
     * @return description
     */
    private static double probD2(double[] cProbOld, int[] k, int i) {

        double result = 0.0;

        if (k[i - 1] == 0) {
            result = 2;
        } else if (((k[i - 1] >= 1) && (k[i - 1] <= 5)) || (k[i - 1] == 10) || (k[i - 1] == 15)
                   || (k[i - 1] == 20)) {
            result = -0.5;
        } else {
            result = 1.0 / 8.0;
        }

        return result;
    }

    /**
     * Method description
     *
     * @param cProbOld description
     * @param k description
     * @param i description
     *
     * @return description
     */
    private static double kProbSingle(double[] cProbOld, int[] k, int i) {

        double product = 1;

        for (int index = 0; index < k.length; index++) {

            if (index != (i - 1)) {
                product *= GHoodHelper.lociProb(k[index], cProbOld[index]);
            } else {
                product = product * probD(cProbOld, k, i);
            }
        }

        return product;
    }
}
